The Challenge of Modelling Local Convective Extreme Events - Does Increasing the Resolution Help?¶
1.0 The study domain and the model domain setup¶
In [7]:
#
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
topocmap = plt.cm.gist_earth
topocmap.set_bad('w', alpha=0)
#First the greater area map
with nc(os.path.join(geodataf,'topo_2.nc')) as topof:
Gtopo=np.ma.masked_equal(topof.variables['topo'][:],0)
Gtopolon=topof.variables['longitude'][:]
Gtopolat=topof.variables['latitude'][:]
#Now only the tiwi-islands manp
with nc(os.path.join(geodataf,'Topo.nc')) as topof:
topo=np.ma.masked_equal(topof.variables['topo'][:],0)
topolon=topof.variables['longitude'][:]
topolat=topof.variables['latitude'][:]
with nc(CPOLF) as fnc:
lon=fnc.variables['longitude'][:]
lat=fnc.variables['latitude'][:]
fig = plt.figure()
ax = fig.add_subplot(111)
ls = LightSource(azdeg=315, altdeg=45)
M = Basemap(projection='rotpole', llcrnrlat=min(Gtopolat), llcrnrlon=min(Gtopolon), urcrnrlat=max(Gtopolat), urcrnrlon=max(Gtopolon),
resolution='c', area_thresh=1, ax=ax, lon_0=0, o_lon_p=0., o_lat_p=90.,)
borders = (([122.963, 139.055], [-17.091, -7.443]),
([127.385, 134.693], [-14.505, -10.005]),
([129.557, 132.529], [-13.757, -10.757]))
text=(('4km res.', 38), ('1.33 km res.', 32), ('444 m res.', 26))
M.pcolormesh(Gtopolon, Gtopolat, ls.hillshade(Gtopo[:], vert_exag=1), cmap= topocmap, ax=ax)
#M.drawmapscale(M.boundarylons[0]-0.7, M.boundarylats[0]+0.7, max(topolon), min(topolat), 10,
# barstyle='fancy', fontsize=19, labelstyle='simple')
for conf in range(3):
lons,lats = borders[conf]
txt, fts = text[conf]
x, y = M([lons[0], lons[0], lons[1],lons[1]], [lats[0], lats[1], lats[1], lats[0]])
xy = list(zip(x,y))
ttt = plt.text(*M(lons[0]+0.1, lats[0]+0.05), txt, fontsize=fts)
ttt.set_bbox(dict(facecolor='white', alpha=0.7, edgecolor='none'))
poly = Polygon(xy, alpha=1, edgecolor='k', lw=3, facecolor='none')
plt.gca().add_patch(poly)
#M.drawcoastlines()
M.scatter(*M(131.043101, -12.250323), marker='*',color='b', s=550)
plt.text(*M(borders[0][0][0]+0.15, borders[0][-1][1]-0.4), 'b)', fontsize=32)
axins = zoomed_inset_axes(ax, 7, loc=2)
axins.set_xlim(min(topolon), max(topolon))
axins.set_ylim(min(topolat), max(topolat))
tmap = Basemap(projection='cyl',llcrnrlat=min(lat), llcrnrlon=min(lon), urcrnrlat=max(lat), urcrnrlon=max(lon),
resolution='c', area_thresh=1, ax=axins)
tmap.pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap= topocmap)
mark_inset(ax, axins, loc1=3, loc2=4, fc="none", ec='k')
plt.annotate('a)', xy=(0.01 ,0.95), xycoords='axes fraction', fontsize=32)
lons, lats = [min(Gtopolon), max(Gtopolon)], [min(Gtopolat), max(Gtopolat)]
axins2 = inset_axes(ax, width='15%', height='15%', loc=3)
Worldmap= Basemap(projection='ortho', lat_0=Gtopolat[int(len(Gtopolat)/2)], lon_0=Gtopolon[int(len(Gtopolon)/2)],
ax=axins2, height=240*10**3)
Worldmap.drawcoastlines(color='white')
Worldmap.fillcontinents(color='gray')
bx, by = Worldmap(M.boundarylons, M.boundarylats)
xy = list(zip(bx,by))
poly = Polygon(xy, alpha=1, edgecolor='none', lw=0, facecolor='darkcyan')
Worldmap.ax.add_patch(poly)
Out[7]:
2.0 Climatology of rainfall and the importance of the diurnal cycle¶
To investigate the links between Monsoon stage and Rainfall extremes in the considered area it is worthwhile having a look at the general climatology of rainfall and investigate the importance of the diurnal cycle in the stduy domain.
The relative importance of the diurnal cycle is investigated in two ways:
- Calculating the fraction of total variability each frequency for the spectral density of the domain avg. rain-rates
- Presenting a map the fraction of variability for the 1/24h frequncy band of the spectral density at each grid-point
- Presenting the diurnal cycle
Calculate the fraction of total spectral variance¶
Because the cpol data contains holes (missing data) we chose to create a Lomb-Scargle periodigram rather than a a fourier decomposition of the data. Lomb-Scargle periodigrams are created by least-spaure fitting sinosoids to the data. The LombScargle Class provided by astropy's stats module will create the periodogram. For best practice we create a function that calculates the fraction of the spectral variance.
In [14]:
#Plot the output map
fig = plt.figure()
#ax1 = fig.add_subplot(111)
ax1 = fig.add_subplot(131)
ax1.set_title('rainfall climatology (1998 - 2017)', size=28)
im = m.pcolormesh(lon, lat, mask*np.nanmean(cpol24h,axis=0)*24, vmin=0, cmap=colmap.GMT_drywet, ax=ax1)
m.drawcoastlines()
cbar = m.colorbar(im,location='bottom',pad="5%")
cbar.set_label('Avg. rain-rate [mm/d]', size=28)
cbar.ax.tick_params(labelsize=24)
#'''
ax2 = fig.add_subplot(132)
ax2.semilogx(f, psd * 9, lw=4)
ax2.set_xlim(0.25,60)
ax2.set_ylim(0,0.6)
ax2.set_xlabel('Periods ($\\tau$) [d]',size=28)
ax2.set_ylabel('Spectral variance [%]',size=28)
ax2.tick_params(labelsize=28)
ax3 = fig.add_subplot(133)
im2 = m.pcolormesh(lon, lat, np.ma.masked_invalid(spec_map)*2 , vmin=0, vmax=0.8,cmap='Blues', ax=ax3)
m.drawcoastlines()
cbar = m.colorbar(im2,location='bottom',pad="5%")
cbar.set_label('Spectral variance of $\\tau$ = 1 d [%]', size=28)
cbar.ax.tick_params(labelsize=24)
# %%
fig.subplots_adjust(right=0.99, bottom=0.83, top=0.99,left=0.01, hspace=0.05, wspace=0.1)
#'''
2.1 The diurnal cycle in different Monsoon stages¶
In [19]:
#Plot the diurnal Cycle
fig, ax = plt.subplots(1, figsize=(12,6))
ax.set_ylabel('Rain-rate mm/d', fontsize=22)
ax.set_xlabel('Local-time', fontsize=22)
ax.plot(idx, dc_burst.values, 'fuchsia', linewidth=2.0, label='Bursts')
ax.plot(idx, dc_break.values, 'g', linewidth=2.0, label='Breaks')
ax.plot(idx, (dc_break+dc_burst)/2, 'firebrick', linewidth=2.0, label='All')
ax.xaxis.set_major_formatter(matplotlib.dates.DateFormatter('%_H'))
ax.legend(loc=0)
ax.grid(linewidth=0)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
2.2 Local Extemes and the stage of the Monsoon¶
The above plots show that diurnal cycle is the most important mode of variability in the considered area. The island thunderstrom hector is the most extreme convective system that occurs at the time-scale at perdiods of 24 h. The Hector 'signal' is most pronounced during the Monsoon break period. Before studying Hector we split the CPOL area into break and burst period and investigate the occurrence of extremes during Monsoon break and burst.
First we define extreme events by applying the 99th percentile threshold of the total 10min dataset. Then we define the break and burst events by the definition of Narsey et al. and calculate the number of burst and the number of break events during the CPOL era
In [21]:
# Now Plot the maps
maxval =[(0,15), (0,15), (-12,12), (0, 150), (0, 150), (-200, 200)]
names=['Avg Rain during breaks', 'Avg Rain during bursts', 'Difference (breaks - bursts)',
'Extremes during breaks', 'Extremes during bursts', 'Difference (breaks - bursts)']
cmap_list = [(colmap.GMT_drywet,'k'), (colmap.GMT_drywet,'k'), (colmap.GMT_polar_r,'k'),
(mpl.cm.magma,'w'), (mpl.cm.magma,'w'), (colmap.GMT_polar_r,'k')]
C = [Mean[1]*1.25, Mean[2], Mean[1]*1.25 - Mean[2], D[1], D[2], D[1] - D[2]]
fig = plt.figure()
mask = ((D[0] * 0 )+1)
im=[]
for i, idx in enumerate(C):
ax = fig.add_subplot(2,3,i+1)
data = np.zeros([len(lat),len(lon)])
colm, coast = cmap_list[i]
try:
im.append(m.pcolormesh(lon, lat, mask*C[i].filled(0),vmin=maxval[i][0],vmax=maxval[i][1],cmap=colm))
except AttributeError:
im.append(m.pcolormesh(lon, lat, mask*C[i],vmin=maxval[i][0],vmax=maxval[i][1], cmap=colm))
m.drawcoastlines(color=coast)
#cbar=m.colorbar(im,location='bottom',pad='5%')
#cbar.set_label(u'Freq. of events above 99th% [%]',size=18)
#cbar.ax.tick_params(labelsize=18)
ax.set_title('%s'%names[i],size=28)
fig.subplots_adjust(right=0.95, bottom=.45, top=0.95,left=0.05, hspace=0, wspace=0)
cbar_ax1 = fig.add_axes([0.05, 0.72,0.58, 0.01])
cbar1 = fig.colorbar(im[0], cax=cbar_ax1, orientation='horizontal')
cbar1.ax.tick_params(labelsize=24)
cbar1.set_label(u'Rain-rate [mm/d]',size=28)
cbar_ax2 = fig.add_axes([0.67, 0.72,0.27, 0.01])
cbar2 = fig.colorbar(im[2], cax=cbar_ax2, orientation='horizontal')
cbar2.ax.tick_params(labelsize=24)
cbar2.set_label(u'Difference',size=28)
cbar_ax3 = fig.add_axes([0.05, 0.47,0.58, 0.01])
cbar3 = fig.colorbar(im[3], cax=cbar_ax3, orientation='horizontal')
cbar3.ax.tick_params(labelsize=24)
cbar3.set_label(u'N. of events above 99th%',size=28)
cbar_ax4 = fig.add_axes([0.67, 0.47,0.27, 0.01])
cbar4 = fig.colorbar(im[-1], cax=cbar_ax4, orientation='horizontal')
cbar4.ax.tick_params(labelsize=24)
cbar4.set_label(u'Difference',size=28)
2.0 Ensemble Simulations of Islands Thuderstorms¶
To investigate how well Extreme Events caused by local tropicsl Island Thunderstorms are represented in a numerical model we choose a period of 8 consecutive occurring Hector events over the Tiwi Island during the Monsoon break period in November 2006. To increase the sample size for more robust statistics we create ensemble of 8 members. The ensemble is created by 8 different initialisation times.
2.1 A Week of Hector¶
Now we compare the ensemble timeseries of the simulated rainfall for the 1.33km and 0.44km Simulations with the CPOL observations.
In [24]:
#Plot area avg TS
import matplotlib.dates as mdates
myFmt = mdates.DateFormatter('%d/%m %H LT')
#
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(20,10), dpi=72)
ax = fig.add_subplot(111)
obs_time_1 = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
obs_time_2 = pd.DatetimeIndex(OBS.dataset[2].coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
um044_time = pd.DatetimeIndex(UM044ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
um133_time = pd.DatetimeIndex(UM133ens.coords['t'].values).tz_localize(utc).tz_convert(timezone).tz_localize(None).to_pydatetime()
## First plot the area avg of the cpol data
#ax.plot(OBS.dataset[0].coords['t'], np.nanmean(OBS.dataset[0].variables['lsrain'].values[:,0,:]*mask.filled(np.nan),axis=(1,2)),
# label='CPOL', color='firebrick')
#ax.plot(obs_time_2, OBS.dataset[2].variables['lsrain'].mean(axis=(1,2)), label='Cmorph 8km',
# color='fuchsia', linestyle='-')
ax.plot(obs_time_1, OBS.dataset[1].variables['lsrain'].mean(axis=(1,2)), linestyle='-',
color='lightseagreen', label='CPOL')
#############
ax.plot(um133_time, UM133ens.variables['lsrain'].mean(axis=(0,2,3,4)), color='fuchsia',
label='UM 1.33km ensmean', linestyle='--')
#ax.fill_between(UM133ens.coords['t'].values,
# UM133ens.variables['lsrain'].mean(axis=(2,3,4)).min(axis=0).values,
# UM133ens.variables['lsrain'].mean(axis=(2,3,4)).max(axis=0).values, color='fuchsia',
# alpha=0.2)
############
ax.plot(um044_time, UM044ens.variables['lsrain'].mean(axis=(0,2,3,4)), color='fuchsia',
label='UM 0.44km ensmean')
ax.fill_between(um044_time,
UM044ens.variables['lsrain'].mean(axis=(2,3,4)).min(axis=0).values,
UM044ens.variables['lsrain'].mean(axis=(2,3,4)).max(axis=0).values, color='cornflowerblue',
alpha=0.2)
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
ax.set_ylim(0.0,0.42)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=18)
ax.xaxis.set_major_formatter(myFmt)
ax.grid(color='r', linestyle='-', linewidth=0)
2.1.1 Maps of Rainfall¶
Now lets look at some animations of UM rainfall occuring over the Tiwi-islands and compare it to the cpol observations
In [54]:
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Ens-2.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
Out[54]:
Maps of Rainfall (0.44km)¶
In [27]:
#Create the avg. maps of rainfall 0.44km
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.35, top=0.95,left=0.05, hspace=0, wspace=0)
UM = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m[0].drawcoastlines()
cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['Obs']
ax[-1].set_title('Obs', fontsize=24)
for i in range(1,9):
ax.append(fig.add_subplot(3,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=100, ax=ax[-1]))
m[-1].drawcoastlines()
#m[-1].shadedrelief()
im.append(m[-1].pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append(datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime('%e %b %H UTC'))
ax[-1].set_title(tit[-1],fontsize=24)
Maps of Rainfall (1.33km)¶
In [28]:
#Create the avg maps of rainfall for 1.33km
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.35, top=0.95,left=0.05, hspace=0, wspace=0)
UM = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01) / 6
um_data = np.ma.masked_less(UM.variables['lsrain'][:].values, -0.01) / 6
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM.coords['lat'][:], UM.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data,vmin=0.0,vmax=5,cmap=colm)]
m[0].drawcoastlines()
cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['Obs']
ax[-1].set_title('Obs', fontsize=24)
for i in range(1,9):
ax.append(fig.add_subplot(3,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=100, ax=ax[-1]))
m[-1].drawcoastlines()
#m[-1].shadedrelief()
im.append(m[-1].pcolormesh(lon2, lat2, um_data[i-1][0],vmin=0.0,vmax=5,cmap=colm))
tit.append(datetime.strptime(ensembles[i-1], '%Y%m%dT%H%MZ').strftime('%e %b %H UTC'))
ax[-1].set_title(tit[-1],fontsize=24)
Ensemble mean and std¶
In [29]:
#Plet the avg maps for comparison
fig = plt.figure()
fig.subplots_adjust(right=0.94, bottom=0.01, top=0.6,left=0.01, hspace=0, wspace=0)
UM_2 = UM133ens.groupby('t.day').sum(dim='t').mean(dim='day')
UM_1 = UM044ens.groupby('t.day').sum(dim='t').mean(dim='day')
obs = OBS.dataset[1].groupby('t.day').sum(dim='t').mean(dim='day')
colm = colmap2.Blues
colm.set_under('w', alpha=0)
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
obs_data = np.ma.masked_less(obs.variables['lsrain'][:].values, -0.01)
data = [obs_data[1:-1]/6,
np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).mean(axis=0)[0]/6,
None,
np.ma.masked_less(UM_1.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0],
np.ma.masked_less(UM_2.variables['lsrain'][:].values/6, -0.01).std(axis=0)[0]
]
data.append(None)
data.append(data[0] - data[1])
data.append(data[0] - data[2])
UM_1min = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_2min = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).min(axis=0)[0]
UM_1max = np.ma.masked_less(UM_1.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
UM_2max = np.ma.masked_less(UM_2.variables['lsrain'][:].values, -0.01).max(axis=0)[0]
#data.append(None)
#data.append(UM_1max - UM_1min)
#data.append(UM_1max - UM_1min)
lat1, lon1 = obs.variables['latitude'][:,0], obs.variables['longitude'][0,:]
lat2, lon2 = UM_1.coords['lat'][:], UM_1.coords['lon'][:]
ax = [fig.add_subplot(3,3,1)]
m = [Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='f',
area_thresh=1, ax=ax[0])]
#m[0].pcolormesh(topolon, topolat, ls.hillshade(topo[:], vert_exag=1), cmap='gray', ax=ax[0])
im = [m[0].pcolormesh(lon1, lat1, obs_data/6,vmin=0.0,vmax=25/6,cmap=colm)]
m[0].drawcoastlines()
#cbar_ax = fig.add_axes([0.14, 0.33, 0.74, 0.01])
#cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
#cbar.ax.tick_params(labelsize=24)
#cbar.set_label('Avg. Rain-rate [mm/d]',size=24)
tit = ['CPOL avg.', 'UM 0.44km ens avg.', 'UM 1.33km ens avg.',
None, 'UM 0.44 ens std.', 'UM 1.33km ens std.',
None, 'CPOL - UM 0.44 ens avg.' ,'CPOL - UM 1.33km ens avg.',
None, 'UM 0.44 ens range', 'UM 1.33km ens range']
colors = {1:colm, 2:colm, 4:colm, 5:colm,
7:colmap.GMT_polar_r, 8:colmap.GMT_polar_r,
10:colm, 11:colm}
Range = {1:(0.0, 25/6), 2:(0.0, 25/6), 4:(0,8/6), 5:(0,8/6), 7:(-20/6,20/6), 8:(-20/6,20/6), 10:(0,30/6), 11:(0,30/6)}
height = {2:0.43, 5:0.23, 8:0.03, 11:0.05}
ax[-1].set_title('CPOL avg.', fontsize=24)
cbar_ax = [fig.add_axes([0.95, height[2], 0.01, 0.15])]
cbar = [fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical')]
cbar[-1].ax.tick_params(labelsize=24)
cbar[-1].set_label('Rain-rate [mm/d]',size=24)
for i in range(1,len(data)):
if data[i] is not None:
ax.append(fig.add_subplot(3,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=100, ax=ax[-1]))
m[-1].drawcoastlines()
#m[-1].shadedrelief()
im.append(m[-1].pcolormesh(lon2, lat2, data[i],vmin=Range[i][0],vmax=Range[i][1],cmap=colors[i]))
ax[-1].set_title(tit[i],fontsize=24)
if i in (5, 8, 11):
cbar_ax.append(fig.add_axes([0.95, height[i], 0.01, 0.15]))
cbar.append(fig.colorbar(im[-1], cax=cbar_ax[-1], orientation='vertical'))
cbar[-1].ax.tick_params(labelsize=24)
cbar[-1].set_label('Rain-rate [mm/d]',size=24)
2.1.2 The Simulated Diurnal Cycle¶
In [25]:
# Plot avg diurnal cycle
#Get the minimum time period that is covered by all simulations and also observations
fig = plt.figure(figsize=(20,10), dpi=72)
ax = fig.add_subplot(111)
####################
cpol_1 = OBS.dataset[1].groupby('t.hour').mean(dim='t').mean(dim=('y','x'))
obs_S = pd.Series(cpol_1.variables['lsrain'].values, index=(cpol_1.variables['hour'].values+9) % 24 ).sort_index(inplace=False)
ax.plot(obs_S.index, obs_S.values, linestyle='-', color='lightseagreen', label='CPOL')
##################
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).mean(dim='ens')
um_2 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).min(dim='ens')
um_3 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).max(dim='ens')
um_044 = pd.DataFrame({'mean':um_1.variables['lsrain'].values,
'min' :um_2.variables['lsrain'].values,
'max' :um_3.variables['lsrain'].values},
index=(um_1.variables['hour'].values + 9) % 24).sort_index(inplace=False)
ax.plot(um_044.index, um_044['mean'].values, color='fuchsia', linestyle='-', label='UM 0.44km ensmean')
###################
um_4 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).mean(dim='ens')
um_5 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).min(dim='ens')
um_6 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('lat','lon','surface')).max(dim='ens')
um_133 = pd.DataFrame({'mean':um_4.variables['lsrain'].values,
'min' :um_5.variables['lsrain'].values,
'max' :um_6.variables['lsrain'].values},
index=(um_4.variables['hour'].values + 9) % 24).sort_index(inplace=False)
ax.plot(um_133['mean'].index, um_133['mean'].values, color='fuchsia', linestyle='--', label='UM 1.33km ensmean')
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
#####################
ax.fill_between(um_044.index, um_044['min'], um_044['max'], color='cornflowerblue', alpha=0.2, label='UM 0.44km ens. spread')
ax.fill_between(um_133.index, um_133['min'], um_133['max'], color='coral', alpha=0.2, label='UM 1.33km ens. spread')
#ax.set_ylim(0.0,0.42)
ax.set_xlim(0,23)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=18)
ax.set_xlabel('Local time [h]', fontsize=18)
ax.legend(loc=0, fontsize=18)
ax.tick_params(labelsize=18)
ax.grid(color='r', linestyle='-', linewidth=0)
In [55]:
#Plot Animatoin
import io
import base64
from IPython.display import HTML
outvid=os.path.join('Vid','WeekOfHector-Diurnal-2.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
Out[55]:
Location and timing of daily rainfall peaks¶
In [48]:
#Plot the maps
plt.close()
um_1 = UM044ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
um_2 = UM133ens.groupby('t.hour').mean(dim='t').mean(dim=('surface'))
obs = OBS.dataset[1].groupby('t.hour').mean(dim='t')
time = (um_1.hour + 9) % 24
data = [obs, um_2, um_1]
#fig = plt.figure(figsize=(25,15))
lat1, lon1 = OBS.dataset[1].variables['latitude'].values[:,0], OBS.dataset[1].variables['longitude'].values[0,:]
lat2, lon2 = um_1.coords['lat'][:], um_1.coords['lon'][:]
#fig.subplots_adjust(right=0.94, bottom=0.025, top=0.67,left=0.01, hspace=0, wspace=0)
first = True
val_max = 3
colm = colmap2.Blues
colm.set_under('w', alpha=0)
m, ax, im = [], [], []
for i in range(24):
break
fname_1 = os.path.join(outdir, 'DiurnalCycle_%02i.png'%i)
if first:
for col in range(3*2):
if col != 3:
ax.append(fig.add_subplot(2,3,col+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1, ax=ax[-1]))
m[-1].drawcoastlines()
###MEAN
ax[0].set_title('CPOL',fontsize=22)
ax[1].set_title('UM 1.33km mean', fontsize=22)
ax[2].set_title('UM 0.44km mean', fontsize=22)
ax[-2].set_title('UM 1.33km std', fontsize=22)
ax[-1].set_title('UM 0.44km std', fontsize=22)
im.append(m[0].pcolormesh(lon1, lat1, obs['lsrain'].values[i], shading='gouraud', vmin=0.1, vmax=val_max, cmap=colm))
im.append(m[1].pcolormesh(lon2, lat2, um_2['lsrain'].mean(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
im.append(m[2].pcolormesh(lon2, lat2, um_1['lsrain'].mean(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
###STD
im.append(m[-2].pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
im.append(m[-1].pcolormesh(lon2, lat2, um_2['lsrain'].std(dim='ens').values[i], vmin=0.1, vmax=val_max, cmap=colm))
first = False
cbar_ax = fig.add_axes([0.14, 0.0, 0.74, 0.02])
cbar=fig.colorbar(im[-1], cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
else:
###MEAN
im[0].set_array(obs['lsrain'].values[i].ravel())
im[1].set_array(um_2['lsrain'].mean(dim='ens').values[i].ravel())
im[2].set_array(um_1['lsrain'].mean(dim='ens').values[i].ravel())
###STD
im[-2].set_array(um_2['lsrain'].std(dim='ens').values[i].ravel())
im[-1].set_array(um_1['lsrain'].std(dim='ens').values[i].ravel())
cbar.set_label('Rain-rate at %2i LT [mm/h]'%((i+9)%24),size=24)
fig.savefig(fname_1, bbox_inches='tight', format='png', dpi=72)
#break
day = (time[1:5] - 9) % 24
night = (time[5:5+4] - 9) % 24
obs_diff = obs['lsrain'][day].sum(dim='hour') - obs['lsrain'][night].sum('hour')
um_1_diff = um_1['lsrain'][day].sum(dim='hour') - um_1['lsrain'][night].sum('hour')
um_2_diff = um_2['lsrain'][day].sum(dim='hour') - um_2['lsrain'][night].sum('hour')
mask = np.ma.masked_invalid(obs['lsrain'].values) * 0 + 1
obs_max = (((np.argmax(obs['lsrain'].values * mask, axis=0)*mask[0] + 9 ) % 24)/2).round(0) * 2
um_1_max = ((um_1['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9) % 24 / 2).round(0) * 2
um_2_max = ((um_2['lsrain'].argmax(dim='hour').mean(dim='ens').round(0) + 9) % 24 / 2).round(0) * 2
colm = CubeYF_20.get_mpl_colormap(N=24, gamma=2.0)
#colm = qualitative.Paired[12].get_mpl_colormap(N=24, gamma=2.0)
m, ax = [], []
data = [obs_max, um_1_max, um_2_max]
names = ['CPOL', 'UM 1.33km ens. mean', 'UM 0.44km ens. mean']
fig = plt.figure(figsize=(25,15))
fig.subplots_adjust(right=0.94, bottom=0.025, top=0.67,left=0.01, hspace=0, wspace=0)
cbar_ax = fig.add_axes([0.12, 0.17, 0.74, 0.02])
for i in range(3):
ax.append(fig.add_subplot(1,3,i+1))
m.append(Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2),
resolution='f', area_thresh=1, ax=ax[-1]))
m[-1].drawcoastlines()
try:
im = m[-1].pcolormesh(lon1, lat1, data[i], vmin=0, vmax=24, cmap=colm)
except TypeError:
im = m[-1].pcolormesh(lon2, lat2, data[i], vmin=0, vmax=24, cmap=colm)
ax[-1].set_title(names[i], fontsize=22)
cbar=fig.colorbar(im, cax=cbar_ax, orientation='horizontal',ticks=list(range(1,24,2)))
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Local Time [h]', size=24)
3.0 Simulated Extreme Events¶
3.1 Pdf's of Rainfall¶
In [37]:
#Plot the pdfs
from mpl_toolkits.axes_grid.inset_locator import inset_axes
dist=lognorm([LogNorm.UM044['shape']*3.4], loc = np.log(LogNorm.UM044['scale']*1.9))#, shape # mu, sigma
med = np.median(um044vec.loc[um044vec>0].values)
X = np.linspace(0,65,200)
#b = 2.5
nbins = 100
def fitpdf(x, b):
a=med*100
return ((b/a) * (x/a)**(b-1)) / (1+(x/a)**b)**2
#popt, pcov = curve_fit(powerlaw, um044cnt['Rain-rate'].values[1:],
# um044cnt['count'].values[1:])
popt, pcov = PowerFit.UM044['popt']/4, PowerFit.UM044['pcov']
fig = plt.figure()
ax = fig.add_subplot(111)
ax2 = ax.twinx()
#Define histogram data to plot
histdata = (um044vec.loc[um044vec>0], um133vec.loc[um133vec>0], obs_vec.loc[obs_vec>0])
facecolors = ['fuchsia', 'firebrick', 'lightseagreen']
#n1_1, bins1_1, patches1_1 = ax2.hist(, bins=nbins, normed=1, facecolor='firebrick', alpha=0.25)
#n2, bins2, patches2 = ax2.hist(, bins=nbins, normed=1, facecolor='lightseagreen', alpha=0.25)
ax2.set_ylim(0,1.4)
ax2.set_yticks([])
ax.plot(np.linspace(0.01, 100,200), fitpdf(np.linspace(0.01,100,200),popt[0]), lw=4, color='navy',
label='Log-logistic $\\beta = %02.2e$'%popt[0])
ax.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
ax.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
ax.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#ax.scatter(um044cnt['Rain-rate'].values, um044cnt['count'].values, c='red', s=125, alpha=0.24)
#ax.scatter(obscnt['Rain-rate'].values, obscnt['count'].values, c='purple', s=125, alpha=0.24)
#ax.plot(np.linspace(0, 50,200), um044pdf, lw=4, label='UM 0.44km ens', color='lightseagreen')
#ax.plot(np.linspace(0, 50,200), obspdf, lw=4, label='CPOL', color='fuchsia')
ax.plot(np.linspace(0.0,50,200), dist.pdf(np.linspace(0.0, 50,200)), color='firebrick', alpha=0.8,
lw=4, label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
n1, bins1, patches1 = ax2.hist(histdata, bins=nbins, color=facecolors, normed=1, alpha=0.25)
#fitted_pl.plot_pdf(ax=ax,color='k',lw=4)
#ax.plot(np.linspace(0.1, 100), um044dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 0.44km area-avg')
#ax.plot(np.linspace(0.1, 100), um133dens_aavg(np.linspace(0.1,100)), lw=4, label='UM 1.33km area-avg')
plt.subplots_adjust(bottom=0.05, right=0.9, top=0.4)
ax.legend(loc=4, fontsize=28)
ax.tick_params(labelsize=28)
ax.set_ylim(0.,.4)
ax.set_xlim(0.,40)
ax.set_ylabel('Density', fontsize=28)
ax.set_xlabel('Rain-rate mm/h', fontsize=28)
axin = inset_axes(ax, width = "60%", height = "60%", loc=1)
axin.plot(cnts.UM044['Rain-rate'].values, cnts.UM044['counts'].values,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
axin.plot(cnts.UM133['Rain-rate'].values, cnts.UM133['counts'].values,'--', label='UM 1.33 km ens',lw=4, color='fuchsia')
axin.plot(cnts.CPOL['Rain-rate'].values, cnts.CPOL['counts'].values,'-', label='CPOL',lw=4, color='lightseagreen')
#axin.plot(X, um044pdf,'-', label='UM 0.44 km ens',lw=4, color='fuchsia')
#axin.plot(X, obspdf,'-', label='CPOL',lw=4, color='lightseagreen')
axin.plot(np.linspace(0.01,80,200), dist.pdf(np.linspace(0.01, 80,200)), lw=4, color='firebrick', alpha=0.8,
label='Log-norm ($\\mu=%02.2f, \\sigma=%02.2f$)'%(scale, shape))
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=24)
axin.set_xlim(0.01, 90)
axin.set_ylim(0.000025, 1)
axin.grid(color='r', linestyle='-', linewidth=0)
axin.annotate("b)", (0.01,0.03), xycoords="axes fraction", fontsize=28)
ax.grid(color='r', linestyle='-', linewidth=0)
axin.plot(np.linspace(0.01, 50,200), fitpdf(np.linspace(0.01,50,200),popt[0]), lw=4, color='navy')
axin.set_ylabel('$\\log$(Density)', fontsize=24)
axin.set_xlabel('$\\log$(Rain-rate)', fontsize=24)
Out[37]:
3.2 Percentiles¶
Now get the percentiles for the ensemble simulation and compare it against the observations
In [39]:
#Plot the percentages
fig=plt.figure()
ax = fig.add_subplot(111)
ax.plot(PERC.index, PERC['Obs'].values, color='lightseagreen', linestyle='-', label='CPOL',lw=3)
ax.plot(PERC.index, PERC['UM 1.33km'].values, color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
ax.plot(PERC.index, PERC['UM 0.44km'].values, color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
ax.set_xlim(1,100)
#ax.set_ylim(10,100)
#ax.set_yscale('log')
#ax.set_xscale('log')
ax.set_xlabel('Percentile [%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)
axin = inset_axes(ax, width = "80%", height = "60%", loc=9)
axin.plot(PERC.index[1:], PERC['Obs'].values[1:], color='lightseagreen', linestyle='-', label='CPOL',lw=3)
axin.plot(PERC.index[1:], PERC['UM 1.33km'].values[1:], color='fuchsia', linestyle='--', label='UM 1.33km',lw=3)
axin.plot(PERC.index[1:], PERC['UM 0.44km'].values[1:], color='fuchsia', linestyle='-', label='UM 0.44km', lw=3)
axin.set_yscale('log')
axin.set_xscale('log')
#axin.set_xticks([])
#axin.set_yticks([])
axin.tick_params(labelsize=24)
axin.set_xlim(50, 100)
axin.set_ylim(.25,100)
#axin.set_xticks([20, 200, 500])
axin.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
#axin.get_xaxis().get_major_formatter().labelOnlyBase = False
axin.grid(color='r', linestyle='-', linewidth=0)
axin.annotate("b)", (0.01,0.03), xycoords="axes fraction", fontsize=28)
axin.set_xticks(range(50,100,10), minor=False)
ax.grid(color='r', linestyle='-', linewidth=0)
axin.set_ylabel('Rain-rate[mm/h]', fontsize=24)
axin.set_xlabel('Percentile [%]', fontsize=24)
Out[39]:
Z. Intermezzo (Comparison CPOL old vs new)¶
In [ ]:
##Intermezzo
first = True
fig = plt.figure()
fig.subplots_adjust(right=0.95, bottom=0.05, top=0.95,left=0.05, hspace=0, wspace=0)
data1 = np.ma.masked_less(OBS.dataset[-1].variables['lsrain'][:].values, 0.02)
data2 = np.ma.masked_less(OBS.dataset[1].variables['lsrain'][:].values, 0.02)
t1 = OBS.dataset[-1].variables['t'][:]
t2 = OBS.dataset[1].variables['t'][:]
tsteps = pd.DatetimeIndex(OBS.dataset[1].coords['t'].values).tz_localize(utc).tz_convert(timezone).to_pydatetime()
#for tt, tstep in enumerate(tsteps):
tt = 0
fname = os.path.join(outdir, 'Comparison_%05i.png'%tt)
#sys.stdout.flush()
#sys.stdout.write('\r Creating %s .... '%os.path.basename(fname))
#sys.stdout.flush()
lon1, lat1 = OBS.dataset[-1].variables['longitude'][0,:], OBS.dataset[-1].variables['latitude'][:,0]
lon2, lat2 = OBS.dataset[1].variables['longitude'][0,:], OBS.dataset[1].variables['latitude'][:,0]
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
m1 = Basemap(llcrnrlat=min(lat1), llcrnrlon=min(lon1), urcrnrlat=max(lat1), urcrnrlon=max(lon1), resolution='i',
area_thresh=1, ax=ax1)
m2 = Basemap(llcrnrlat=min(lat2), llcrnrlon=min(lon2), urcrnrlat=max(lat2), urcrnrlon=max(lon2), resolution='i',
area_thresh=1, ax=ax2)
im1 = m1.pcolormesh(lon1, lat1, data1[tt].filled(-1),vmin=0.0,vmax=30,cmap='Blues')
im2 = m2.pcolormesh(lon2, lat2, data2[tt].filled(-1),vmin=0.0,vmax=30,cmap='Blues')
m1.drawcoastlines()
m2.drawcoastlines()
cbar_ax = fig.add_axes([0.14, 0.35, 0.74, 0.01])
cbar=fig.colorbar(im1, cax=cbar_ax, orientation='horizontal')
cbar.ax.tick_params(labelsize=24)
cbar.set_label('Rain-rate [mm/h]',size=24)
ax2.set_title('Old CPOL version at %s LT' %tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
ax1.set_title('New CPOL version at %s LT'%tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
def ani(tt):
im1.set_array((data1[tt]).filled(-1).ravel())
im2.set_array((data2[tt]).filled(-1).ravel())
ax2.set_title('Old CPOL version at %s LT' %tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
ax1.set_title('New CPOL version at %s LT'%tsteps[tt].strftime('%d. %b %Y %H:%M'), fontsize=24)
return(im1, im2)
# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, ani, frames=100, interval=20, blit=True)
#fig.savefig(fname, bbox_inches='tight', format='png', dpi=72)
#break
#sys.stdout.write('\n')
In [ ]:
#Visualize intermezzo
import io
import base64
from IPython.display import HTML
viddir = os.path.join(os.getenv('HOME'),'Data/Extremes/CPOL/Plot/Versions')
outvid=os.path.join(viddir,'WeekOfHector.mp4')
video = io.open(outvid, 'r+b').read()
encoded = base64.b64encode(video)
HTML(data='''<video alt="test" width="950" height="500" controls>
<source src="data:video/mp4;base64,{0}" type="video/mp4" />
</video>'''.format(encoded.decode('ascii')))
2.1 Distribution of Percentiles during break and bursts¶
Now we compare the distributions of percentiles during break and burst periods. We also add
In [ ]:
#Plot percentiles for breaks, bursts and all
colors = {'10min':'firebrick'}# '1h':'lightseagreen', '3h':'fuchsia' , '6h':'cornflowerblue'}
fig=plt.figure()
ax = fig.add_subplot(111)
for con in (('break','-'), ('burst','--'), ('all','-.')):
for gr in colors.keys():
df = pd.read_hdf(PercF, con[0])[gr]
ax.plot(df.index, df.values, color=colors[gr], linestyle=con[-1], label='%s (%s)'%(gr, con[0]))
ax.set_xlim(1,99.9999)
#ax.set_ylim(0.11,)
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_xlabel('Percentile [%]', fontsize=26)
ax.set_ylabel('Rain-rate [mm/h]', fontsize=26)
ax.legend(loc=0, fontsize=24)
ax.tick_params(labelsize=24)